% conv.M
% compute recursive posterior mean 
% based on Metropolis-Hastings output

% Taeyoung Doh
% created: 08/29/2008
%--------------------------------------------------------------------1.0
% house cleaning
%--------------------------------------------------------------------
clear all;
close all;
diary off;

%--------------------------------------------------------------------
% default parameters
%--------------------------------------------------------------------

% call environments and model selection
rsmodel_spec;
%model_specsv; 
loaddata;

%--------------------------------------------------------------------
% prior specification
%--------------------------------------------------------------------
global prior_;
prior_ = feval(fnames_.prior,msel);

nobs    = size(data,1); 

global datasel psel 
% MHRUN
%  length of mhrun defines the number of columns in graph table
%	odd integer  (e.g., 1) for constant regime 
%	even integer (e.g., 2) for regime shifts 
 


block1   = 151;		% starting block
blockn   = 250;		% ending block
nsim     = 10000;	% number of simulation draws in each block

hpdprob  = 0.90;	% probability for highest posterior density

% for data density (modified harmonic mean)
hmax    = 20;

%--------------------------------------------------------------------
% model descriptions and basic calculations
%--------------------------------------------------------------------

% parameter names
pnames = feval(fnames_.pnames,msel);

%--------------------------------------------------------------------
% means and standard deviations
%--------------------------------------------------------------------

mhrun = 1;
mhrunid  = num2str(mhrun','%2.0f');
 

% load path
 %lpath = [path_.para 'mhrunpmchain2' mhrunid  '/' ];
  lpath = [path_.para 'mhrunpm' mhrunid '/' ];  
% initialize output
ndraws     = 0;
drawmean   = 0;
drawsqmean = 0;
sumdrawsq  = 0;
mpos       = zeros(blockn-block1+1,1); 
rho_inf    = zeros(2000,2); 
drawblock  = []; 
t=1; 

% loop for block
for jblock=block1:blockn

	% load file
	lfile = [lpath 'mhdraw' num2str(jblock,'%04.0f')];
	load(lfile);

	% use only p.d. covariance matrix results for prior draws
	if mhrun==0
		parasim = parasim(find(retcode==0),:);
	end
     
	% number of simulations in each block
	nsimul = size(parasim,1);

	% collect simulations
    
    drawblock = [drawblock;parasim]; 
    drawdim   = size(drawblock,2);

	% compute sums of x(i) and x(i)^2 to be used to calculate means and s.d.
	ndraws     = ndraws + nsimul;
end

% compute means and standard deviations
drawmean   = cummean(drawblock);

% save output
% sfile = [path_.para 'conv_mhrun' mhrunid '.mat'];
% save(sfile,'drawmean');

% plot posterior expected probability of regime 

%mm = mod(drawdim,4); 
%if mm == 0;
%   parax = 4*(drawdim/4-1)+1;
%else 
%   parax = 4*floor(drawdim/4)+1; 
%end    

 %  for i=1:4:parax
 %      figure;
 %      j=1; 
 %      while j<5;
 %          if i+j-1<=drawdim
 %      subplot(2,2,j); plot(drawblock(:,i+j-1));
 %          j=j+1;
 %          else
 %          j=5;
 %          end 
 %       end 
 %  end 

figure 
drawblock = drawblock(1:100:end,:); 


for j=1:16 
    xxn  = deblank(pnames(j,:));
   subplot(4,4,j); plot(drawblock(:,j)); xlabel(xxn,'FontSize',8);
end 
figure 
for j=17:drawdim 
   
    xxn  = deblank(pnames(j,:));
     subplot(4,4,j-16); plot(drawblock(:,j)); xlabel(xxn,'FontSize',8);
end 

% Split sample mean test 

nmp = 0.5*length(drawblock);
p   = 2; 
thetam = zeros(drawdim,p);
 
spec = zeros(drawdim,1); 
tst  = zeros(drawdim,1); 

for j=1:drawdim 
           
       
          lag = 0.1*nmp;
       %  lag =0.01*nmp; 
    
        thetam(j,1) = mean(drawblock(1:nmp,j)); 
        thetam(j,2) = mean(drawblock(nmp+1:end,j)); 
         factor1 = var(drawblock(1:nmp,j));
         factor2 = var(drawblock(nmp+1:end,j)); 
         for jj=1:lag
covx1     = cov(drawblock(1:nmp-jj,j),drawblock(1+jj:nmp,j));
covx1     = covx1(1,2);
covx2     = cov(drawblock(nmp+1:2*nmp-jj,j),drawblock(nmp+1+jj:2*nmp,j));
covx2     = covx2(1,2);
factor1  = factor1  + 2*(lag+1-jj)*covx1/(lag+1);
factor2  = factor2  + 2*(lag+1-jj)*covx2/(lag+1);
         end    
        spec(j) = factor1 + factor2;    
   
    tst(j) = sqrt(nmp)*(thetam(j,2)-thetam(j,1))/(sqrt(spec(j)));
  
end

 tst 

